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Our method of computing lattice sums 1 is carried through the final stage 
of numerical compulation. We calculate the expansion coefficients of the 
crystal potential, evaluated at each inequivalent site of a number of cubic 
structures. Some of the structures chosen are simple, to afford comparison 
with previous calculations. Many, however, are too complex to be amenable 
to other methods, and our results are the first reported. The sequence of 
computer programs, timing, and accuracy are discussed. 

In a previous paper 1 concerning lattice sums we discussed summation 
methods which hinge on two facts: First, many lattices can be decom- 
posed into so-called "primitive" lattices. Second, a wide class of summa- 
tion procedures is feasible for such primitive lattices, which are im- 
possible or impractical for arbitrary lattices. In particular, we developed 
an extension of the Ewald method, following the general philosophy of 
Nijboer and DeWette, 2 Adler, 3 and Barlow and Macdonald. 4 The condi- 
tions under which a given lattice can be decomposed into primitive lat- 
tices, and the algorithm accomplishing this, have been discussed with 
full generality and rigor by Graham. 5 

In this paper we present a sample of numerical results obtained by 
our method, together with some discussion of computational techniques 
and computational efficiency. We have evaluated the coefficients, up to 
order six, for the spherical harmonic expansion of the potential due to a 
lattice of point charges. We have done this for every inequivalent site 
in a number of cubic lattices. Some of the lattices were chosen to pro- 
vide comparison with other calculations; some, on the contrary, 
were chosen because then- complexity puts them beyond the reach of 
other methods; others still because cubic lattices without an inversion 
center are interesting in a number of other connections. 
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We define terms and display the relevant equations, as briefly as pos- 
sible. For a fuller exposition we refer to Refs. 1 and 5. 

A lattice is primitive if it consists of a set of points of position r„ and 
charge q n such that 



r„ = 2 riibi, n, = 0, ±1, ±2, 



(la) 



?» = 9o(-D ni+n2+n3 . 



(lb) 



For instance, NaCl is primitive, with bi = (1,0,0), b 2 = (0,|,0), and b 3 
(0,0,5), where the components of the b's are given along the cube axes 
in units of the cube dimensions. We call the vectors b,- the basis vectors 
of the primitive lattice. These should be distinguished from primitive 
translations c, , which are defined by the relation q(i) — g(r + WiCi + 
n£t + W3C3), for all r, and for all integers n\ , n% , n 3 . 

An arbitrary three-dimensional lattice can be decomposed if it is of 
periodicity 2 N X 2 N X 2 N . This condition can be loosened to periodicities 
2 L X 2 M X 2 N with L,M < N, since by taking appropriate multiples 
in the direction in which the periodicity is 2 L or 2 M one can reduce this 
case to the previous one. The existence condition can be further loosened 
to include lattices composed of "interlocking" sublattices of periodicity 
2 L X 2 M X 2^, provided each such sublattice is separately electrically 
neutral. 

We consider, then, an array of periodicity 2 N X 2 N X 2 N , and primi- 
tive translations Ci , C2 , C3 . This array can be decomposed into 2 _ 
primitive lattices. These component lattices will be grouped into 3N 
sets. All lattices within one set will be defined by the same basis vectors, 
but will differ in choice of origin. With respect to the origin in the 
original lattice, the position of the origin of a component lattice is de- 
noted by R. We define auxiliary vectors a,- = Ci/2 N . We wish to express 
the basis vectors bj of the primitive component lattices, as well as the 
origin positions R, in terms of the a,- . Let the columns of the matrix A 
denote the components of a, , and likewise for B and b t - . Superscripts 
will denote sets of primitive component lattices, subscripts denote lat- 
tices within a given set. Then the required relations are 
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(2b) 



.1 1 -1. 

B(n+3) — 2B (n) 

|R<") = (0,0,0) 
{R(«+D| = JR(«)} + {RCt + B^j 

9 („+d = qMy/2. 



(2c) 

(2d) 
(3a) 
(3b) 
(4a) 
(4b) 



These definitions differ from those given in Refs. 1 and 5 by relabeling 
some of the indices, which has the effect that all the b vectors now de- 
fine right-handed bases. The primitive lattices so defined are orthonor- 
mal in the sense that 



Z q a U) (r n )q (i) (T n ) = 8««. 



(5) 



where n runs over the sites in a unit cell, the superscript i or j denotes a 
set of primitive lattices, and the subscript a or denotes a particular 
lattice belonging to that set. The primitive lattices do not necessarily 
correspond to physically real structures. If we take A = 1, representa- 
tive lattices belonging to the first three sets are shown in Fig. 1. There 
is one member in set 1, two in set 2, four in set 3. For every set we show 
only the lattice with origin at R = 0. Only set 3 has a physical counter- 
part, and will be recognized as the CsCl structure. Further examples 
may be found in Ref. 5. 

We now consider the expansion coefficients of the crystal potential. 
Just as the physical lattice is represented as the sum of primitive lat- 
tices, so a given expansion coefficient is computed as the sum of the 
corresponding coefficients proper to the primitive lattices. As suggested 
in the previous paper, 1 the coefficients arising from primitive com- 
ponent lattices can be computed once and for all, so that the problem 
for a real physical lattice, no matter how complicated, is reduced to the 
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Fig. 1 — Conventional unit cell representation of primitive lattices, sets 1, 2, 
and 3 (Translation R = 0,0,0). 

decomposition of that lattice into primitive components. If the potential 
is given, as usual, by 

F(x) = E E Ci m I x \ l Y,,^(e x tVs ) (6) 

1=0 m=-l 

then the expansion coefficients may be represented as 

c- E E c lM (b M) - R« M) ) (7) 

sets (i) translations (a) 

where the notation is consistent with that of (2) through (5). Shortening 
the argument of the primitive lattice Ci m to p, (p = b (i) — R« (i) )> we 
have, for each primitive lattice, 

4tt 



Clm = 



E Qn | Qn | ^.m(^p n) ¥'e n ) 



(8) 



2Z + 1 

where n runs over the sites of that particular lattice. The spherical 
harmonics are defined as usual: 



Yi.m(0#) 



-[■ 



21 + 1 (2- |ml) n 



!]•' 



'P, m (0). 



(9) 



4ir (Z+|m|) 

For small J, the summation in (8) converges poorly or not at all. By 
means of the formalism discussed in Ref. 1, (8) can be transformed into 
the following, which converges rapidly for all I: 

4tt 1 



C ta = 



21 + 1 T(l + |) 

•{Elf.f*r(i + Wli.f) 

n 

•F im (^ n ,^J(-l) ni+n2+ " 3 

+ tW E exp (Ow.-*/ ) | n„ T 2 



(10) 



■r(i,™r"| n„ \*)-Yi m (e nH ,<pn n ) - ^1- 



COMPUTATION OF LATTICE SUMS 2437 

Here, V is the incomplete Gamma function, and v c is the cell volume 
bi • b 2 X b 3 . The vector n in reciprocal space is defined by 

n„ = h„"> - i(h/'» + h,« + h,«) (11) 

where the h's are basis vectors in the reciprocal lattice, 

h« = i (b 2 ,, ' ) X h. w ), (12) 

Vc 

and h„ is a lattice vector of the reciprocal lattice. The sums on n are to 
be understood as triple sums on ?h ,n,2,n 3 . The scale factor a is given by 

a = »r* (13) 

and makes the rate of convergence equal in both coordinate and recip- 
rocal space, regardless of choice of units or size of cell. The result of the 
summation is independent of the choice of a, but the correct choice is 
crucial in practice since it can affect computation time by orders of 
magnitude. The last term in (10) exhibits explicitly the correction dis- 
cussed in Ref . 1 for the case p = 0. 

In practice, there are four computational processes, performing the 
following functions : (i) We generate all the primitive component lattices 
of periodicity 2 N X 2 N X 2 N . We did this for N = 1,2,3, requiring re- 
spectively 7, 63, and 511 lattices, (ii) For each component lattice, we 
calculate the C tm , both real and imaginary parts, with ^ I ^ G, 
<> m ^ I. (Hi) We decompose a given lattice, by taking appropriate 
dot-products (sec (5)). (iv) We reconstitute the Ci m 's for the given 
lattice by combining the results of steps (ii) and (Hi). Parts (i) and (ii) 
are done once only for a given type of geometrical grid. Parts (Hi) and 
(iv) are done for each charge configuration that can be placed on that 
grid. 

The computing time is spent almost entirely on Part (ii). About an 
hour on the IBM 7094 was required to calculate all the Ci m 's for all the 
primitive lattices. To compensate, the application of these results to 32 
physical configurations, for all 56 coefficients, took about twelve seconds. 

The computation for Part (ii) is done shell-wise. By the kth shell we 
mean here a set of points r„ = «ibi -f- n 2 bo -f- w 3 b 3 such that at least one 
rii is equal to k in absolute value, but no n,- is greater than k in absolute 
value. The number of points in the A:th shell is (24fc 2 + 2), and the 
number of points in all shells up to and including the kth. is (2k + l) 3 . 
These shells differ from "Evjen" shells, which, unlike ours, use frac- 
tional charges, and preserve crystal symmetry and electrical neutrality 
for every shell. The feth shell defined above contains an excess charge of 
2<?( — l)*i if the charge q( — 1)" 1+ " 2+ " 3 is associated with each site. Since 
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the excess charge alternates in sign, and the fractional excess decreases 
as k 2 , convergence is not affected adversely. Summation was stopped 
when the contribution of a shell was less than 10 -6 of the aggregate sum. 
This required, on the average, about 4 shells. Allowing for normal round- 
off errors, we expect an accuracy, in final results, of five significant fig- 
ures. 

The computed expansion coefficients are exhibited in Tables I through 
XIII. These tables are reproductions of computer output. Each table 
is headed by the chemical formula, and by the length of the unit cell, 
in angstroms. Next follows a list of all atomic positions in the unit cell. 
The crystallographic data is taken from Wyckoff. 6 Next follows a tabular 
listing of the expansion coefficients, i.e., the C tm of (6). Due to evolution 
in notation, these coefficients are called ALM in Tables I-XIII. Both 
the real and imaginary parts are given for positive m; our definition of 
the spherical harmonics requires C*,_„, = C*, m *. The coefficients are 
listed for the expansion of the potential at every inequivalent site. For 
purposes of this computation, sites whose environment differs only by 
rotations and/or reflections are considered equivalent. All the structures 
are cubic, except SnF 4 . This substance in tetragonal, with c = 7.93 A, 
a = 4.048 A. Since c is so nearly twice a, we have cheated a little by 
stacking together four of the tetragonal cells to make an almost cubic 

Table I 

NA-CL A=5.640 
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Table II 



CL 0.500 0.5.00 0.500 
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5 
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6 
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Al M REAL 


ALM [MAG 


-1.7*998032 


0. 
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Table III 

A=5.*09 
0.500 0.500 0.500 0. 0. 0.500 0.500 

SO 0.750 0.750 0.250 0.750 0.750. 0.750 0.250 
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Table IV 

CAF2 A=5.463 

0- 0. 0. 0. 0.500 0.500 0.500 0. 0.500 0.500 0.500 0. 
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0.250 
0.250 



0.250 
0.250 



0.250 
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4 
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0. 
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0. 
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a. 
o. 
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0. 
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0. 
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0. 
0. 



0. 

0. 

0. 

0. 
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0. 

0. 

0. 



0. 
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Table V 



CU2-B A=4.270 
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0. 
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0.01274866 


0. 
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Table VI 
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Table VII 
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Table VIII 



AC 

a 



0.250 0.250 0.250 



AG2-03 A=4.940 
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cell of average dimension 8 A. We have included SnF 4 mainly as an 
illustration of a case where the actual periodicity is not the same in all 
dimensions. 

For purposes of checking, there exist simple relations between the 
Modelling constant and our coefficients Coo ■ The Madelung constants 
for NaCl, CsCl, ZnS, CaF 2 , Cu 2 hark back to an extensive tabulation 
by J. Sherman. 7 (Sherman's number for Cu 2 0, quoted throughout the 
literature, e.g., by Born and Huang, 8 is in error. The correct value has 
been given by Hund. 9 Very accurate computations for some cubic crys- 
tals have been made by Benson and Zeggeren. 10 For MX compounds, 
the Madelung constant is related to our C 00 by a simple multiplicative 
factor. This factor is (4*-)* times the ionic charge divided by the length 
of the unit cell. For compounds like CaF 2 and Cu 2 0, the Madelung con- 
stant is related in the same fashion to the average of the absolute values 
of the Coo's for the various sites. For Cu 2 0, the potential at each site 
separately has more recently been computed by Dahl. 11 Madelung cal- 
culations for the perovskites have been considered by Templeton, 12 
Fumi and Tosi, 13 and Cowley. 14 Cowley has calculated numerically the 
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"partial M addling coefficients" as defined by Templeton, and Fnmi and 
Tosi have shown how the potential at any site of an arbitrary perovskite 
structure can be obtained by linear combinations of these partial Made- 
lung coefficients. Our numerical comparison for the perovskite Coo's is 
therefore, essentially with Cowley's calculation. Comparison of our Coo's 



Table XIV — Comparison between Present Computation and 
Previous Work for the Coefficients C o 
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Present Result 


Previous Result 


Source for Previous Result 


NaCl 


either 


2.19679 


2.1968 
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either 


1 . 74998 


1.75003 


Benson, Zeggeren 


ZnS 


either 


4.95817 


4.95845 


Benson, Zeggeren 


CaF, 


1 Ca 1 + | F | 


7.55090 


7.55090 


Benson, Zeggeren 


Cu 2 


| Cu | + | 1 
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Hund 


Cu 2 


Cu 


3.14084 
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Dahl 
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5.3768 
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2.9765 


Cowley 
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6.1711 


Cowley 
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Sr 
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4.8905 


Cowley 


SrTi() 3 


Ti 
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with previous work is made in Table XIV. In general, some divergence 
appears in the fifth significant figure. 

We have ourselves checked a portion of our results by using the Evjen 
summation method. This method is particularly useful for coefficients 
with higher /-values, but converges quite slowly for Coo . Only for CsCl 
and NaCl were we able, in reasonable computing time, to carry the 
summation to a sufficient degree of convergence for C 00 . The Evjen 
results confirm our Ewald results for these substances, for all the coef- 
ficients. For the perovskites, we pushed the summation at the Ta site 
as far as the ninth Evjen shell (24,564 neighbors), but for Coo shell-to- 
shell fluctuations were still in the fourth significant figure. For higher 
/-values, the Evjen summation was able to confirm our Ewald results 
for all coefficients, for all sites, for both perovskites. The same is true 
for Cu 2 0. It should be pointed out that the computation of the coeffi- 
cients (exclusive of Coo) takes about ten minutes of 7094 time per crystal 
if the Evjen method is used. We attempted a check on one site of 
Al 2 MgO,i . The computer was stopped after 20 minutes, at which time 
usefully convergent results were obtained only for coefficients with 
/ ^ 3. 

We note that the appearance of coefficients of odd I, at any site, is 
associated with lack of inversion symmetry. 

We believe that the computations we have discussed show that in 
certain applications our method of doing crystal sums offers definite 
advantages, both in accuracy and in time. Suitable applications are 
those where a number of structures must be considered that can be de- 
scribed by basically the same coordinate grid, especially if the number 
of ions per unit cell is large. 
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